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Many-body instabilities of the half-filled honeycomb bilayer are studied using weak coupling renor- 
malization group as well as strong coupling expansion. For spinless fermions and assuming parabolic 
degeneracy, there are 4-independent four-fermion contact couplings. While the dominant instabil- 
ity depends on the microscopic values of the couplings, the broken symmetry state is typically a 
gapped insulator with either broken inversion symmetry or broken time reversal symmetry, with a 
quantized anomalous Hall effect. Under certain conditions, the dominant instability may appear in 
the particle-particle (pairing) channel. For some non-generic fine-tuned initial conditions, weak cou- 
pling RG trajectories flow into the non-interacting fixed point, although generally we find runaway 
flows which we associate with ordering tendencies. Additionally, a tight binding model with nearest 
neighbor hopping and nearest neighbor repulsion is studied in weak and strong couplings and in each 
regime a gapped phase with inversion symmetry breaking is found. In the strong coupling limit, 
the ground state wavefunction is constructed for vanishing in-plane hopping but finite inter-plane 
hopping, which explicitly displays the broken inversion symmetry and a finite difference between the 
number of particles on the two layers. Finally, we discuss the spin-1/2 case and use Fierz identities 
to show that the number of independent 4-fermion contact couplings is 9. The corresponding RG 
equations in the spin-1/2 case are also presented, and used to show that, just as in strong coupling, 
the most dominant weak coupling instability of the repulsive Hubbard model (at half-filling) is an 
anti-ferromagnet . 



I. INTRODUCTION 

The problem of interacting fermions on the A — B 
stacked honeycomb bilayer at half-filling has attracted 
attention due to a confluence of several factors. First, 
purely on theoretical grounds, in its simplest form with 
the nearest neighbor hoping only, the tight-binding ap- 
proximation gives rise to a band structure with two bands 
touching quadratically at the Fermi levels near two non- 
equivalent points in the Brillouin zone, K and K'. Even 
at the non-interacting level, such quadratic degeneracy 
gives rise to logarithmically divergent susceptibilities 2,3 
in several channels as temperature, or frequency, are 
taken to zeroir— . As a result, some form of spon- 
taneous symmetry breaking is expected at finite tem- 
perature upon inclusion of even weak interactions 2 - - — . 
And while fine-tuning is necessary to achieve such band- 
structure, in that (with the exception of square checker- 
board and Kagome lattices studied in Ref4) inclusion 
of trigonal warping termsi eventually gives rise to four 
Dirac fermions at each K-point, non- interacting suscep- 
tibilities may be sufficiently enhanced that many-body 
instabilities appear, albeit at finite coupling strength. 
In this sense, the A-B stacked honeycomb bilayer prob- 
lem is another example of the observation that there 
are no generic weak coupling particle- hole instabilities 9 . 
Rather, fine-tuning, in the form of nesting for example, 
is necessary to bring the strong coupling physics down 
to weak coupling. If we are interested in accessing the 
symmetry-broking phases in the particle-hole channel, as 
we are in this case, then fine tuning is a small price to pay 
for this access, made available within perturbative RG. 
Second, the isolation of graphene bilayers and the exper- 



imental ability to perform, for example, electrical^— , 
angle resolved photoemission^ 3 ., Raman spectroscopy^ 
or infra-red^ measurements, while controlling the gate 
voltage through the neutrality point, gives rise to the 
opportunity to test such theoretical expectations in a 
reasonably well controlled physical setting. In addition, 
the technological promise of this material fuels further 
need to understand its electronic structure and with it 
the many-body interactions. Finally, the problem of in- 
teracting fermions on the AB-stacked honeycomb bilayer 
may soon be realized in cold atom optical lattices, where 
the theory may also be tested. 

The issue of band-structure fine-tuning notwithstand- 
ing, the type of leading instability in a graphene bilayer 
(with spin 1/2 fermions) has been a subject of debate 
as well. A mean-field approach has been used to argue 
for an insulating state with broken inversion symmetry 7 . 
A similar approach has also been argued to lead to triv- 
ial gapped insulating phases^ as well as to an anomalous 
quantum Hall phase 16 . On the other hand, the lead- 
ing weak coupling instability can be analyzed without 
resorting to uncontrolled approximations by using weak 
coupling renormalization group. This approach was used 
in Ref4 where a nematic phase was found to be the dom- 
inant instability within the model studied. Such insta- 
bility was subsequently also argued for in Ref. 6 -. On the 
other hand, an inversion symmetry breaking insulating 
phase has been claimed in Ref A 

To determine what type of broken symmetry state is 
preferred in the case of spinless fermions, we perform 
weak coupling RG analysis by studying the flow of 4 inde- 
pendent symmetry allowed short-range interactions. We 
find that generically, depending on the initial values of 



the 4-fermion contact couplings, the system flows into 
a gapped phase with either broken inversion symmetry 
and a finite difference between the total number of parti- 
cles on the two layers, or broken time reversal symmetry. 
The former state was not found to be preferred in the 
model for spin-1/2 fermions studied in Ref.— (where the 
nematic state was found to dominate), but an example 
of the latter state corresponded to one of the fixed points 
found therein. In particular, for the spinless case studied 
here, we find that a gapped state with anomalous (zero 

2 

B-ficld) quantum Hall conductivity ±2^- has the most 
divergent susceptibility for a range of initial couplings 
as determined by the (right) sink of the RG trajectories 
shown in Fig. ([3]). While non-generic, we also specify spe- 
cial conditions under which the interacting model flows 
back to the non- interacting fixed point. 

In addition, we analyze the specific microscopic model 
with nearest neighbor hopping(s) t (and t±) and nearest 
neighbor repulsion V in both the weak coupling RG and 
in strong coupling. In both regimes we find the (trivial) 
insulating phase with broken inversion symmetry to dom- 
inate. As discussed in more detail below, in weak cou- 
pling the RG flow tends to the left sink shown in Fig. ([3]), 
with a susceptibility that dominates over other broken 
symmetry states mainly due to subdominant terms. In 
strong coupling, we construct a ground state wavefunc- 
tion for V > 0, V± > 0, t = 0, but t± ^ 0, which shows 
explicitly the broken layer inversion. Since in this model, 
the same symmetry appears to be broken in the limit of 
both weak and strong coupling, it is reasonable to as- 
sume that such a broken symmetry state appears at any 
V, V± > 0. 

A similar analysis is presented in the spin-1/2 case 
with short range interactions. For the repulsive Hubbard 
model, we find that the most dominant weak coupling 
instability is towards an anti-ferromagnetic state. Since 
the same ordering tendency happens in the strong cou- 
pling, it is reasonable to assume that in this model, the 
Neel ordering appears at any U > 0. 

This paper is organized as follows: in Section II we 
write down the (non-interacting) bilayer Hamiltonian 
first in the tight-binding approximation and then within 
k • p perturbation theory. In Section III we construct 
the low energy effective theory at the neutrality point 
by fine-tuning the trigonal warping terms to zero. The 
rest of that section deals with identifying microscopic- 
symmetry-allowed 4-fermion contact interaction terms 
using the method of Herbut, Juricic and Roy— used for 
the same purpose in single-layer graphene. Before the re- 
duction due to Fierz identities, there are 9 such couplings 
which further reduce to 4 once Fierz identities are taken 
into account. The weak coupling RG is presented in Sec- 
tion IV, along with the flow diagram in the space of cou- 
pling constant ratios and the analysis of the susceptibility 
growth. The t—V model with weak and strong coupling 
limits is studied in Section V. In Section VI, the spin- 
1/2 case is revisited. Symmetry is used to construct an 
eighteen-dimensional Fierz vector along with the 18 x 18 



Fierz matrix to show that there are 9 independent cou- 
plings in this case. Their RG equations are determined 
and while more general, they are shown to reduce to the 
ones studied in Ref£ under conditions outlined therein. 
In Section VII we study the Hubbard model in weak and 
strong coupling. Section VIII is devoted to conclusions. 
Details of the derivation are presented in the Appendices. 



II. BILAYER HAMILTONIAN 

In this section we will define the non-interacting model 
by using two different approximation methods. First, the 
well known tight binding approximation- will be used 
and then the k- p-method, or equivalently the method of 
invariant s 1 ' 6 ' 18 ! 19 . Both methods lead to the same form 
of the low energy Hamiltonian and it is ultimately a ques- 
tion of convenience which one should be adopted. 




(i) (it) 



FIG. 1: (i) Schematic representation of the A-B stacked bi- 
layer. The low energy wavefunction near K is also sketched, 
with w — e l7V ^ z = i + i 2 ^- The primitive lattice vectors are 
Ri = \/3ax and R2 = ^-ax + \ay. The area of the unit cell 

is A uc — z- (Ri x R 2 ) = a 2 , (ii) Schematic representation 
of the (reciprocal) k-space. 



A. Tight-binding approximation 

The non-interacting Hamiltonian in the tight-binding 
approximation can be written as 

H = H± + H% (1) 
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The low energy field theory is written in terms of the 
eight-component Fermi fields (two layers, 1 and 2, two 
valleys, K and — K, and two sublattices a and b as 
sketched in Fig.([T])): 
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+ ^W^W + u^(r)^i(r) ) . (11) 



In the case of bilayer graphene, the values of the hopping 
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integrals t, t±, t ± were extracted experimentally in 
Ref. 20 . If we define the Fourier transform of a Fermi field 
as Cj{r) — N U c^ 2 J2k e * rc j,k where c = a or b, j = 1 or 
2, and N uc is the number of unit cells. Next, we let 
X t = (a\ k , a\ k , b\ k , k ^ to write the non-interacting 
Hamiltonian <jTJ) as 
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In the above, the wavevector dependent function d^ = 
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\ya, and £3 = — ya. Near K, dic+k 



Y^g e tk ' S where the sum runs over Si = ^-xa 
82 = —^-xa 



-£t{kx + iky) = -v F k + . Near -K, <2_K+k ~ ^tftx ~ 
iky) — vpk-. The low energy spectrum of this (well- 
known) Hamiltonianii2i, which is easily diagonalized, 
will be discussed in the next section. 



B. k ■ p approach 

Instead of resorting to the tight-binding approxima- 
tion, we can also arrive at the low energy Hamilto- 
nian by analyzing the symmetry of the bilayer poten- 
tial alone. This is a standard technique when dealing 
with semiconductors^ and one which has also been ap- 
plied to graphene^. For the sake of self-inclusiveness, we 
present this method as well to show that one arrives at 
the same general form of the Hamiltonian as in the tight- 
binding approximation, although in practice the coeffi- 
cients of various symmetry-allowed terms must be deter- 
mined from experiment. We start with the Schrodinger 
equation for a particle moving in potential due to the 
atoms in layers 1 and 2 separated by 2c 
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The rapidly-varying Bloch functions at K and at K' = 
— K are related by complex conjugation, uk(i") = 
u*_ K {r), irrespective of the layer or sublattice index. 

Moreover, the Bloch functions u^(r) and u^(r) trans- 
form irreducibly under point group operations of the lat- 
tice (see FigfT]). For the sake of concreteness, within the 
nearly free electron approximation for electron wavefunc- 
tions |xi,2) confined to layers 1 and 2 respectively we have 
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i.e. the interlayer hopping arises from the mixing of the 
sublattices a\ and a 2 . The matrix elements of the in- 
plane momentum operator p are also dictated by sym- 
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gives us the effective Hamiltonian near K to read 
t± v 2 k- \ 
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This is equivalent to what we found in the tight-binding 
approximation. 

The spectra of the k • p and the tight-binding Hamilto- 
nians are well known and have been discussed extensively 
in the literature (See e.g~&22-). In the vicinity of each 
K-point, there are four Dirac points: one isotropic at ±K 
and three anisotropic ones arranged in accordance with 3- 
fold lattice symmetry around the isotropic one. When we 
neglect trigonal warping terms, by setting v± = Vi = 0, 
or set the higher order hopping terms t± = = 0, the 
four Dirac points merge into a parabolic degeneracy. 



III. LOW ENERGY EFFECTIVE THEORY 

In the weak coupling limit, the kinetic energy dictates 
which modes are important to determine the behavior of 
the system at low energies. Clearly, at k = we have two 
degenerate levels and two levels at ±ij_. Since we wish 
to work with a theory for the low energy modes only, we 
need to project out the bands which originate from the 
two "split-off' bands. We can do so in several equivalent 
ways. The method used here implements the path inte- 
gral formalism, where we integrate out the Fermi fields 
associated with a\ and a 2 modes (sites), and arrive at 
an effective action with an effective "Hamiltonian" for 
the low energy modes. In addition to the wave vector 
dependence, this "Hamiltonian" is frequency dependent 
as well. Near the K-point, the effective quadratic action 
after integrating out the a-modes is 

e - S eff = e ~ fi dr^ b '[d^+H bb ]^ b 
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Since the integral is Gaussian, we can easily perform it 
and find that up to an additive constant Sjtyf = 
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Within the k • p theory, the parameters V\ and v 2 should 
be determined from experiment. To make contact with 
the notation in literature, Refill have v% = v% and v 2 — 
— V4 (see their Eqs. 6 and 15). 



If we are interested in the modes near the Fermi level 
of an unbiased bilayer, we can simply set uj n — in the 
effective action (|I9j|22[) . As will be obvious from the dis- 
cussion in the next section, terms arising from the cor- 
rections are perturbatively irrelevant near the Gaussian 
fixed point in the sense discussed in a different context 
in Ref^. 

In what follows we will also set v% — v 2 = to fine- 
tune the system to quadratic degeneracy. Such a situ- 
ation arises if in the tight-binding formulation we con- 
sider only the nearest neighbor hopping integrals, t and 
tj_. Otherwise, as mentioned in the introduction, the ul- 
timate low energy dispersion involves four (one isotropic 
and three anisotropic) Dirac conesi&2I. While such fine- 
tuning appears artificial, it is an example of the maxim 9 
that there are no generic weak coupling particle-hole in- 
stabilities. Rather, fine-tuning, in the form of nesting 
for example, is necessary to bring the strong coupling 
physics down to weak coupling. If we are interested in 
accessing the symmetry-breaking phases in the particle- 
hole channel, as we are in this case, then fine tuning is a 
small price to pay for this access made available within 
perturbative RG 9 . 

Putting back the — K point, the low energy degrees of 
freedom can now be expressed in terms of a four compo- 
nent Fermi field 

tft (r) = ^ Jt (r),^ t (r),^ t (p),^ t (r)) , 

i.e. the electronic degrees of freedom are expanded as 
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The non-interacting low energy (imaginary time r) La- 
grangian, which includes both K and K' valleys, and 
which will serve as our (gaussian) fixed point of depar- 
ture, can therefore be written as 

Co = f d 2 r V f (r,r) (|-+ £ Z a d a p ) ^(r, r)(24) 

where we defined the vector function dk and the 4x4 
matrices YI X ' V as 
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The effective mass parameter entering the above equa- 
tions is m — t±/(2v F ) (In the tight binding approxima- 
tion v F — it/ {2a)). The four component Fermi objects 
tj) appearing in Eq. (1241) were defined as the envelope 
Fermi fields in Eq. (|23|) . In the above, the first Pauli ma- 
trix acts in the valley ±K-space and the second in the 
layer 1,2-space. To make contact with the literature we 



also use Dirac 7-matrices which we represent as 
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The action J drCo is invariant under the scale trans- 
formation 
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This means that the " dynamical critical exponent" z = 2 
for the gaussian theory, which will be our point of depar- 
ture when analyzing weak coupling instabilities. 



A. Short range interactions 

From the above discussion of the gaussian fixed point, 
it is evident that the short range interactions, when pro- 
jected onto our low energy modes, will contain among 
other (perturbatively irrelevant) terms, contact four- 
fcrmion terms which are marginal by power counting. 
The rest of this section deals with identifying such 
symmetry-allowed interaction terms. The method used 
here follows almost verbatim the method used by Herbut, 
Juricic and Roy^i in their analysis of the short range in- 
teractions in single layer graphene. In addition to the lat- 
tice symmetries used in Ref.— , we also include the three- 
fold rotational symmetry^, which reduces the number of 
independent four-fermion couplings to 4. 

We can therefore start by writing the general La- 
grangian 



where Cq was introduced in Eq. (|24p and 
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where, at this point, the sum over S includes all six- 
teen independent four-by four matrices (generators of 
SU(A)) and so does the sum over T. Naively, we have 
16 + 8 * 15 = 136 couplings to consider. Just as in the 
case of the single-layer graphene^, this number will be 
dramatically reduced first by using the discrete symme- 
tries of the lattice and second by using Fierz identities. 

The key role in this reduction is played by the be- 
havior of the Bloch functions u(r) under symmetry op- 
erations, which dictates the transformation properties of 
the four component, slowly varying, envelope Fermi fields 
q/>(r) ls i 19 . The dimer centered rotation by 27r/3, mirror 



reflection about the yz-plane and about the xz-axis fol- 
lowed by xy-plane respectively give 

C 3 tp(x,y) = -e-^^Tp^-^x-^y^x-^ 

= - e %™^^~x-^y,^x-~y^7) 

o-v^(x,y) = nl 2 ip(-x,y) = iju 5 tp(-x,y) (38) 
ala*ij)(x,y) = ba'iV'O) ~v) = l2ip(x, -y) (39) 

The time reversal symmetry and translational symmetry 
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In the above, R = toRi + nR,2 where Ri = \/3ax and 



R 2 = ^ai + ±ay. And since K = ^%-x, K R = 

2 ^-(2m + n), where m,n = 0,±1,±2, .... The lattice 
translational symmetry therefore corresponds to the Z 3 
discrete analog of the chiral £4(1) generated by 7375. 



1. Symmetry reduction 

Following Herbut et.aliU, we split the sixteen linearly 
independent four-by-four matrices S and T into four sets 

A = {l4,72,«7o73,«7i75} (42) 

B = {i7o7i>-«7375,no75,ni73} (43) 

C = {7o ,«7o72, 73^7273} (44) 

D = {71^7172,75^7275}- (45) 

The matrices which belong to the set A are even under 
both reflection operations (|38|) and (|39|) . The matrices in 
the set B are odd under y-reflections |38|) and even under 
"x" -reflections (13T)|) . The matrices in the set C are even 
under y-reflections (|38l) and odd under "x" -reflections 
(f3T)f . And finally, matrices belonging to the set D are odd 
under both (f3"5| and (|3T)1) . This means that only quartic 
terms combining matrices from the same set are allowed 
by symmetry. Each such set contains 4 + 2 * 3 = 10 such 
terms and that leaves 40 couplings. 

Eight matrices A\ 2, B\ 2, C1.2 and Di 2 are left invari- 
ant under the spatial translation operation (|4Tj) . These 
give rise to 3 * 4 = 12 couplings, eight direct gXjXj 
(X = A, B,C, 01D and j — l,or2), as well as four 
mixed gx t x 2 ■ I n addition, there are four sets of pairs 
which transform as vectors under (|4Tj) : a — {^3,^3}, 
(3 = {B 4 ,^ 4 }, 7 = {^3,^3} and 6 = {C A ,D A }. These 
give rise to additional 6 couplings. Schematically, four of 
them are Sp= a ,/9,7,<5 12j=i 9pPj ® Pj and two mixed ones 
are g a p(A 3 ®A i -B 3 ® B 4 ) and g 1 s(C 3 ® C 4 + D 3 ® D 4 ). 
Altogether, after inclusion of the translation symmetry, 
we are left with 18 couplings. 



The starting point is the SU (4) algebraic identity (see 
Eq.(A4) of RefiiS) 



s l3 T mn = -Tr[sr a 7T 6 ]ric 
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which leads to 



FIG. 2: Diagrams appearing at 1-loop RG. Each vertex can be 
represented by either a 4x4 matrix (spinless) or 8x8 matrix 
(spin- 1/2) case. 



The unitary part of the time reversal operations 0, 
Eq. (|40[) , happens to correspond to the mirror reflection 
about y, (Eq 155)) . which has already been taken into ac- 
count. However, complex conjugation, further restricts 
the number of couplings. Specifically, mixed terms with 
one purely real and one purely imaginary matrix cannot 
appear, therefore gci_c 2 = 9DiD 2 = 9-yS — 0. This leaves 
15 couplings^ 7 -. 

The lattice symmetries considered by Herbut et.a£l 
did not contain site- or plaquette- centered rotation^ by 
120°. As stated in Eg. (1571) . this symmetry is generated by 
*7i72- Including this symmetry requires that the cross- 
terms gAiA 2 — 9-Bi-B 2 = 9af) = 0. Moreover, it requires 
that 3^2,42 = 9d 1 d 1 , 9BiBi = 9c 2 c 2 , an d 9p = gs- 

This leaves us with the following 9 terms 



(^(x)S^(x)) (^(y)T^(y)) 
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Tr[Sr a TT b ] (^ (x) rV(y)) 0/%)rV(x))(.48) 



The minus sign comes from ip and ip' being anti- 
commuting (four component) Grassman fields. For con- 
tact terms x = y and the above equation (|48|) constitutes 
a set of linear relations between different terms of our 
symmetry reduced interaction Lagrangian (j46j). 
If we arrange the quartic terms into a vector 

V = {{^A^f,{^A 2 ^f + {^D^) 2 , 

{^Bxi>f + (^C 2 ^) 2 , (^B 2 t(;) 2 , (^dV) 2 , 
{^D 2 ^) 2 1 {^A^f + (^B 3 ipf, 
{^B^) 2 + {^A A ^) 2 + {^d^pf + {^D^f, 
(^C 3 ^) 2 + (^D 3 ^) 2 }, (49) 



(50) 



then the Fierz identities lead to the linear constraint 

g AlAl {^ A^) 2 +g A2A2 [(^A 2 i/j) 2 + D^) 2 ] 
9B lBl [{^B^f + ^C 2 ^f] + 9 b 2 b 2 {^ B 2 ^) 2 
9c c (tp^Ciip) 2 + qd d (ift D 2 ip) 2 A straightforward, though somewhat laborious, applica- 

v, if . ,so , / itn ,\2i , \rii/~i r\2 / ,t n ^9ition of (1481) leads to the explicit form of the Fierz matrix 
g a [(^AsiP) + {^B^) ] + 9l [(^C 3 V0 + ftW) \ n the c ^ of spinless fermions 

gp [{^B^f + {^A^f + {^C^f + {^D^) 2 } .(46) 



The 9 terms can be further reduced to 4 independent 
ones by using Fierz identities. 



2. Fierz identities 

We set gxx — gx to continue with the notation of 
Ref.— . We use the method employed therein to write 
down Fierz identitie a 17 ' 23 which, due to the Grassman 
nature of the Fermi fields, relate various seemingly unre- 
lated couplings. 
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.(51) 



The matrix F has four zero eigenvalues and as a result 
there are four independent couplings 17 . 



In order to make a connection with the previous work 3 , we choose to eliminate 



{^B^) 2 + (^C 2 ^) 2 = 

{^B 2 ^) 2 = 
(^C^j) 2 = 
{^A^) 2 + {^B^) 2 = 



-2{^A^f + [(^A^) 2 + (^D^) 2 ] - 2(^D 2 ^) 2 - 2 [{^C^f + (^D^) 2 ] 

(52) 

- [(^A 2 ip) 2 + {^D^) 2 } + (^D 2 i/j) 2 + [{^C^f + (^D 3 ^) 2 ] (53) 
-(^A^) 2 + [{^A 2 i>) 2 + (^D^) 2 ] - 2(^D 2 i>) 2 - [(V f C 3 V) 2 + (^D 3 ^) 2 ] (54) 



-2{^A^) 2 - 2{^D 2 ^) 2 - [(^C 3 ^) 2 + {^D 3 ^f 



(55) 



and 



-2 [(^A 2 ^f + {^D^f} + 4(V t £' 2 V) 2 
+2 [(^C 3 V) 2 + (V f W) 2 ] , 



(56) 



in favor of the remaining four terms. These equations will be used in deriving our RG equations, since elimination of 
fast modes will generate terms such as, for example, (ip 1 * B 2 ip) 2 ■ The above equations show that such a term does not 
correspond to a new coupling in a renormalized action, but rather is a linear combination of terms already present. 
Finally, we arrive at our interaction Lagrangian 



d 2 r 



g Al AxiP(r,r)) 2 + g A2 ((V' 1 4^0, r)) 2 + (^D 1 i/;(t, t)) 
9p 2 (^D 2 i>(r,T)) 2 +g y ((^C 3 i,(r,T)) 2 +(^D 3 i>(r, T )) 2 



(57) 



Above is the most general four-fermion contact interaction Lagrangian for spinless fermions allowed by the symmetry 
of the A-B stacked honeycomb bilayer. In the next section, we study the weak coupling RG flow of the four couplings 
9 An 9a 2 i 9p 2 an d 9-y- The first three couplings appeared in our previous work 3 - where we called them gi, g 2 , and g 3 . 
The fourth coupling, <j 7 , did not appear there since the starting point assumed only finite 171 and, as we will see later, 
<7 7 is not generated if its starting value is zero. 



IV. RENORMALIZATION GROUP ANALYSIS 



four coupling constants to be 



Clearly, gsr's are marginal by power-counting and the 
question is how they flow. The RG procedure employed 
here follows Refj22 and consists of integrating out the 
fermionic modes in a thin shell between the initial cutoff 
A and A/s, while the integral over u> extends from —00 
to 00. Since we are working in weak coupling, we can 
integrate out the fast modes perturbatively in g's. The 
diagrams needed are shown in Fig.@. Afterwards, the 
lengths r, times r and the modes ip are rescaled accord- 
ing to Eqs. ([32l and the change of the coupling constants 
is noted. (To the order we are working, the dynami- 
cal critical exponent z remains 2). While the details of 
the derivation are provided in the Appendix, we note in 
passing that the analysis is facilitated by the use of the 
identities 



dg Al 
dins 

dg A2 
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dgp 2 
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- / ^9A 1 gA 2 



171 

An 



(60) 



- (9M ~ 2g Al g A2 + 8g\ 2 - 2g A2 g D2 + g 2 D , 2 + 

4 (3gA 2 - gp 2 ) g 1 + 6g 2 ) ^- (61) 

2 (2ffi 2 + 2g Al g D2 - &g A2 gp 2 - ^9p 2 + 
8<m 2 s 7 + 2g 2 ) £ (62) 

771 

-2 9l (~2g Al + 2g A2 + 2 9l ) — . (63) 



i"" dw 
_™ 2ir 



d 2 k 



Gk(iuj) ® G T k(=Rw) 



±u ® u + x 7q ® 7q ) 

V 0=1 / 



771 
47T 



In ; 



where the non-interacting Green's function is 



These equations reduce to the ones studied in Ref.— when 
we set g 1 = in this work and N — 2 in Eqs.(6-8) of 
RefA Their analysis proceeds along the lines discussed 
(58) in Refi 3 . We note that each RG equation corresponds to 
a quadratic polynomial in coupling constants. Therefore, 
dividing each equation by g A2 (which is g 3 in the notation 
of Ref£), we obtain three equations 



Gk(icj) — (—icu + S ■ dk)~ 



(59) 



and, just as before, d£ = d%_ = S K = 72, 

and Yi y =71. 

Using this procedure, we find the RG equations for the 



dgA t 
dgA 2 
dgp 2 
dgA 2 
dg 1 
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9 A t gp 2 g 1 
9A 2 ' 9A 2 ' gA 2 

9 At gp 2 g 1 
9a 2 ' 9a 2 ' 9a 2 
9 At gp 2 g 1 

9A 2 ' 9A 2 ' gA 2 



(64) 
(65) 
(66) 




FIG. 3: Flow diagram in the coupling constant ratio space assuming that gA 2 < (generic behavior). There are two sinks given 
by Eqs. (|73[) and (|75|l in the text and two mixed fixed ratios (|74[1 and ([76}. For <ja 2 > the flow is reversed, and generically, 
the divergent coupling constant ratios simply mean that gA 2 has shrunk and crossed 0. After this point qa 2 becomes negative 
and the directionality shown here is restored. 



where ratios 



j9A 1 



gjj^- = -2*+* u (a* £**l) (to) 

agA 2 9a 2 \gA 2 9A 2 9A 2 J 

TZ-lofx V z) = — - rl 9 ° 2 / \ 

x 2 -2x + 8-2y + y 2 +4{3-y)z + Qz 2 gA 9a 2 _ 9d 2 | ^ / 9A 1 9d 2 9 1 \ ^ 

(67) 2 d 9A 2 9A 2 32 \gA 2 ' 9a 2 ' 9A 2 ) 

„ , , 4 (1 + xy - 3y - y 2 + 4z + z 2 ) d jt 2 9i . ^ ( 9a 1 9d 2 9i \ ,- oA 

fc32{X,y,Z) = \ ; ' ffA 2 "l = 1" #42 , , • (72) 

x 2 - 2x + 8 - 2y + y 2 +4(3 - y) z + 6z 2 dg A2 9A 2 \gA 2 9a 2 9a 2 ) 

Note that the right hand side of these equations is a func- 

2z (— 2x + 2 + 2z) ^j on f coupling constant ratios only, i.e. it is autonomous 

- 2x + 8 - 2y + y 2 + 4 (3 - y) z + 6z 2 ' in the new variables S±, iSz. Jhu, We can think of 

QA 2 QA 2 9A 2 

(69) the right hand sides effectively as (highly non-linear) /3- 
functions for the ratios. The advantage of rewriting the 
flow equations this way, is that in this form it is easier to 
analyze the qualitative nature of the flow diagram. Un- 
Equations (I64H66[) are homogeneous, which means that like in the case where g 7 was assumed to vanish from the 
we can instead study the flow of the coupling constant start 3 , in the present case the /3-function for gA 2 is not 



72.42 (x, y, z) = 



negative semidefinite. It may appear therefore, that we 
lose the directionality of the flow equations in the three 
dimensional ratio space. This turns out not to be the 
case, since the (ellipsoidal) region in the 3D ratio space 
where dgA 2 /d\iis changes sign is precisely the same re- 
gion where the "/3" -functions for the ratios change sign, 
and so it is enough to determine the directionality of the 
flow of the trajectories near fixed points of the ratios, 
which turns out to be simple enough. 

The qualitative analysis proceeds by finding the fixed 
points in the ratio space. There are four of them: 



9 Ax 9p 2 g 7 

9*A 2 9*A 2 9*A 2 

9m 9*d 2 9*~f 

9a 2 9a 2 9a 2 

9\ 9*p 2 9^ 

9*A 2 9*A 2 9*A 2 

9*Ax 9*p 2 9* 

9a 2 9a 2 9a 2 



(0,-1.085,0) 
(0,0.566,0) 
(0,6.519,0) 
(-1,-1,-1) 



(73) 
(74) 
(75) 
(76) 



The first three are the N = 2 analog of the (N — 4) fixed 
ratios found in Ref4, while the fourth one is new. For 
gA 2 < the stability analysis gives the first (|73|) and the 
third one (|75|) to be sinks (see FigJ3]). The second one 
([74)1 is mixed, with two stable directions (negative eigen- 
values) and one unstable direction (positive eigenvalue). 
The fourth one (|76|) is also mixed, with one positive, one 
negative and one zero eigenvalue. For qa 2 > the di- 
rectionality of the flows is reversed and the sinks become 
sources while the mixed fixed points remain mixed but 
also with reversed sense of flow. These "runaway" flows 
in the coupling constant ratio space for gA 2 > simply 
correspond to decrease of gA 2 which eventually crosses 
zero, where the ratios become infinite, and then become 
negative. Once negative, the flows are described by the 
two stable sinks, separated by a critical plane (3D version 
of the (red) separatrix shown in Fig. 3 of Ref£). 

The generic flow for initial gA 2 of any sign is towards 
large and negative gA 2 and towards either one of the two 
ratio sinks. 

However, it is interesting to ask, under which (non- 
trivial) conditions, may all the coupling constants flow 
to zero. One possibility, is to fine-tune the initial values 
of gA x and g 7 to zero, set the initial value of gA 2 > and 
the ratio —1.085 < go 2 lgA 2 < 6.519. In this case, the 
flow is towards both gA 2 — > and gr> 2 — > while their 
ratio approaches 0.566. Note that in this case we have 
to fine-tune two of the four symmetry allowed couplings 
gAiiOi to vanish. 

Another possibility involves the new fixed point in the 
ratio space at ( -^-) = (-1, -1, -1). While the 

1 \ 9A 2 9A 2 9A 2 J ^ ' ' ' 

fixed point is mixed, in that one of the RG eigenvalues 
is negative and one positive, one eigenvalue, whose right 
eigenvector is (-^, vanishes. This means that in 

the vicinity of this fixed point, the flow along or against 
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A A 


0,0,0,0 


1,-1,-4,-2 


0,0,0,0 


1,-1,0,0 


A s 


1,-1,0,2 


0,0,0,0 


0,0,0,0 


1,-1,0,0 


Ac 


2,2,-4, -4 


1,-1,0,2 


2,2,-4,-8 


1,-1,0,0 


A D 


1,-1,-4,-2 


1,-6,-4,4 


2,2,-4,-8 


1,-1,0,0 



TABLE I: The susceptibility coefficients Ax 3 ■, Bx 3 n C Xj n v x j 
in Eg. ([81)) for different particle-hole order parameters ip'Oiip. 



this direction in the ratio space is very slow. Impor- 
tantly, our numerical integration finds that for gA 2 > 0, 
starting anywhere(l) along the line [ ££l £Ea. ] = 

b y w ° \9A 2 1 9A 2 ' 9A 2 ) 

(—A, —A, —A) for < A < 1, the flow is towards A = 1 
with decreasing gA 2 — >• 0. The flow trajectories passing 
through this line segment in the ratio space, however, 
are not straight lines. In fact, they connect with the 

fixed point ( = (0, -1.085,0). This means 

1 \9A 2 ' 9A 2 ' 9A 2 ) v ' ' ' 

that there is a nontrivial (curved) finite surface in the 
ratio space along which the flow is directed towards the 
non-interacting fixed point if gA 2 starts out positive. In 
this case only one parameter needs to be fine-tuned in 
order to start on this surface. This interesting behavior, 
however, is non-generic, in that such a surface is unsta- 
ble, and the generic flow for initial gA 2 > is towards 
large and negative gA 2 and towards the two ratio sinks. 



A. Susceptibilities and ordered states 



The physics associated with the fixed ratios analyzed 
in the previous section can be understood by studying 
the flow of the susceptibilities toward forming various 
orders. For translationally invariant order parameters, 
the susceptibilities can be calculated from the above flows 
by introducing source terms into the action so that S — > 
S + AS 



AS = -A°* J dTd 2 r^(r,T)0^ip(r,T) 

- A° p ; f drd 2 rMr,r)O^Mr,r) (77) 



Next, we integrate out the fermionic modes within a small 
shell given byA/s<fc<A and find the correction to the 
source term perturbatively in the g's. We then substitute 
the flow of the g's into the pref actors of various source 
terms and ask which diverges the fastest as s increases. 



Particle-hole channels: 



In the particle-hole channel, we therefore find: 

A°-( S )40^<(r,r) = .s 2 A^(l)40^< 



S 2 A -(l)^ SMM n A4^<(r,T) 

M 

s 2 A°-(l) Yl .9ma^< T om^< (r, r) (78) 



M 



where i is summed over the 16 independent order param- 
eters (generators of S'C/(4)), and 
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OM 



- OM 



did 


r A d 2 k 


>. 


/a/ s (27t) 2 


duj 


r A d 2 k 




/a/ s (2^) 2 



r B:[G k (ia;)O i G u (w)Ml79) 
MG k (icj)OiG k (ia;)M (80) 



Using (|55|) one can easily convince oneself that the only 
non-zero contributions to Hom come from O — M, and 
that the matrix Tqm is proportional to O. From here 
we find the flow equations for the source terms 

din Ax . , . m 

± = 2+(A Xj g Al + B X]9D2 + C Xj g A2 + V X]9l ) — 

(81) 

where X = A,B,C,D and j = 1,2,3,4. The results 
of this calculation, i.e. the values of A Xj , B Xj , C Xj , and 
T>x - , are shown in Table HI The coupling constants g 
are functions of s and, in order to determine the most 
likely ordering tendency, it is necessary to find out which 
source term Ax,- grows the fastest. We can write each of 
these equations as 



d\nA Xj 
din s 



= 2 



C 



gA 2 



As • Lis 9d > 



X, 



V 



x, 



gA 2 

g 7 \ m 



gA 2 



gA 2 J 4?r 



(82) 



and near the two sinks, we can take gA 2 < and substi- 
tute the fixed point ratios. 

Near the first sink f 4s ^l 2 -, -4-) = (0,-1.085,0) 

V 9 A2 ' 9a 2 9a 2 ) v ' ' > 

and, plugging in these values, we find that the fastest 
divergence appears for AS, A^, and A^ 3 . Discrimi- 
nating between the first and the last two ordering ten- 
dencies requires knowledge of the sign of the subleading 
ratio -^ 2 -. 

9A 2 

If j^- < 0, i.e. if it approaches from below as s 

increases, then the most dominant particle-hole order- 
ing tendency is towards a finite expectation value of 
C\ = 70 = I20V Physically, this order parameter, which 
corresponds to an imbalance in the number of particles 
on the two different layers, opens up a gap at the K and 
— K points in the Brillouin zone and the system is a (triv- 
ial) insulator. As shown below, this turns out to be the 



case for a lattice model with a nearest neighbor repulsion 
V. 

On the other hand, if -^ L - > and approaches zero 

9A 2 

from above as s increases, then the most dominant or- 
dering tendency among the particle-hole channels stud- 
ied here is towards finite expectation values of C3 = 73 = 
r x a y and -D3 = 75 = r v a v , both of which are odd under 
time reversal symmetry (|40[) . 

Near the second sink [ 4^, ^l 2 -, ■$-) = (0,6.519,0) 

\9a 2 ' 9a 2 ' 9a 2 ) _ v ' ' ; 

the most dominant ordering tendency is towards a fi- 
nite expectation value of D2 = 47172 = t z o~ z . The 
corresponding order parameter also opens up a gap in 
the single particle spectrum, but unlike C\, it breaks 
time reversal symmetry. This results in an anomalous 

quantum Hall state, with zero B-field Hall conductivity 

2 

a xy — ±2^-. Such a state is a bilayer analog of the Hal- 
dane model for the quantum Hall effect without Landau 
levels 2 ^. 



2. Particle-particle channels: 

Since our Fermions are spinless, if the integral 

d 2 rMr,r)0 { 2Mr,r) 

is to be finite, we must have O l ab = —0\ a . Of the six- 
teen SU(4) generators (|42|) this condition selects the six 
matrices B3, B4, G2, G3, D\, and D±. Integrating out the 
fast modes, we are left with the following renormalization 
of the source term for the slow modes tp < 

A^( S )^<(r,r)0i^<(r,r) = S 2 A^(l^<0^<(r, 

O l a pG-k t a (ioj n )M ab il) b< (r, r)G_k,«c(M^„)M crf ^ rf< (r,*r 

(i 

Evaluating the necessary matrix products leads to 
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d In s d In s 



m 
-in 
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m 
4^ ( 



(84) 
-(85) 
(86) 
(87) 



Near the first sink, = (0, -1.085, 0), and, 

V 9a 2 9a 2 9a 2 J 

substituting these values into Eqs. (154"l) - (|57| . we find that 
the strongest divergence appears for A^ p 4 . The leading 
divergence is as fast as for Ag, Ag, and A^f, but it 
differs in the subleading terms. In fact, for gA 1 > 0, the 



strongest divergence is in the particle-hole channel 
discussed above. In principle, fine-tuning and keeping 
the subleading term gAi/gA 2 < (and setting g 1 = or 
keeping g n j gA 2 < 0) may lead to the strongest divergence 
appearing in the particle-particle channel. 



V. t - V MODEL FOR SPINLESS FERMIONS 

While the weak coupling results are quite general, we 
can apply them to a specific microscopic model, which 
happens to be quite revealing in that we can also analyze 
it for strong coupling and thus compare the two regimes. 
We consider spinless fermions hopping on the half-filled 
A-B stacked honeycomb bilayer, with nearest neighbor 
hoppings t and t± only and with nearest neighbor repul- 
sions V and V±. The corresponding Hamiltonian is 



2t±/(9t 2 ). The interaction term can be written as 



t~>i,n.t — 



£[(*- 

q.G 



(r) - <2 q rii, 2 ,- q (r)) x 



(^q+G^&x.q+GC?") ~ ^G-q^.G+qC?"))] 



(94) 



a reciprocal lattice vector and rib . = 

T\T~ lA'.k+q- I n addition, each fermionic mode is 

restricted to reside in the first Brillouin zone. 



where G is 

Ek & J,lA,k+q; 

Taking k 

in the above sum to be near K or — K gives two possibili- 
ties for q: either q ~ or q ~ ±2K. Note that in the first 
case do = 3 while in the second case (i± 2 K = 0. There- 
fore, only the first term contributes a marginal coupling, 
and the above Hamiltonian gives rise to the low energy 
interaction Lagrangian 



n = 



hI + v 1 - 



V" (88) 
where 

H o = -t±j2( a i^ a ^ + h - 

R 

Hi = -t^(6 t 1 (R + (5)ai(R)+4(R-^)a2(R) + /i.c 



where 



(89) 



rt-V _ (0) 



(0) 

9oi 



8t x uc ' 



(95) 



(96) 



R,<5 



V 1 - = Vi^^(R) ai (R)-i) (^(R),m ( R) 

v" = y^f4(R) ai (R)-i 

+ F^f4(R)a 2 (R)-i 



The area of the unit cell A uc = z ■ (Ri x R 2 ) = ^^a 2 . 
This means that we should start our RG flow with a small 
and negative (attractive) g^(ip^Ciip) 2 , which should be 
rewritten using the Fierz identity (|54[) . The initial con- 
ditions are therefore 



R,<5 



6 t 1 (R + 5)6i(R + <5) 
bl(R-5)b 2 (R.-S) 



g Al (8 = l) = -g%, g D2 (s = l) = -^ i c ) 1 (97) 
g Aa (s = i) = g<g, g y (s = i) = -gW. (98) 

— Substituting these as the initial conditions into our RG 
^equations we find that none of the coupling constants 
^(change sign and they all diverge at the same value of 
s. The ratios of the couplings flow to the fixed point 

( 4s 4 s -, S-) = (0, -1.085,0). Therefore, as discussed 

in the previous section, the fastest divergence appears in 
the channel l2<r z . We therefore conclude that the weak 
coupling instability of this model is towards a gapped, 
broken inversion symmetry state with an imbalance of 
the number of particles on layer 1 compared to layer 2. 

B. Strong coupling limit 



A. Weak coupling 

In order to project onto the low energy modes, we first 
rewrite the Hamiltonian (155)) as an imaginary time Grass- 
man path integral. We then integrate out the a\ and a 2 
modes perturbatively. This results in 



C 



(t-v) 
eff 



tu- 



rn 



where 



Co = 



C 



ill I 




(93) 

Fourier transforming the Fermi modes in the first term 
gives rise to the kinetic energy term (|24p with m = 



Setting t = t± = we find three ground states at half- 
, \ filling: 

S ,t) + h.c.j (92) ^ each site of sublattice a\ and of 6 2 is singly occupied 

ii) each site of sublattice a 2 and of b\ is singly occupied 
\ 2 iii) each site of sublattice b\ and of 6 2 is singly occupied. 

Each of these states breaks sublattice symmetry, but 
the average density of particles on each layer is the same 
and equal to 1 per unit cell (which contains two sites in 
each layer). The states i) and ii) differ from the state 

iii) by the occupation of the ai-a 2 dimer which is singly 
occupied for the former and empty for the latter. 



n b2 (R - 5, t) 




& 2 |0} 



□ 



□ (cos0al +sin6»a^ |0) 



FIG. 4: Schematic representation of the strong coupling state 
for the spinless t— V model. For t = but finite t± , V , and V± , 
the sublattice &i is empty, while sublattice &2 is fully occupied. 
The oi — a 2 dimer is singly occupied, with an electron partially 
delocalized onto ai despite the repulsion from the occupied 62 
sites. Such gapped state breaks inversion symmetry between 
the layers 1 and 2. 



If we now set t = but t± 7^ then we can further 
lower the energy of i) and ii) by delocalizing the electron 
on the dimer. So, consider the deformation of the state 
i): we seek a state of the form 

= n( cos6lat i( R ) +sin6>4(R))4(R + o)|0;(99) 

R 

For t± — 0, we have = 0, but once t± 7^ we expect 
0^0. 

Acting on \^fg) with H (for t = 0) and requiring \^fg) 
to be an eigenstate gives: 



V± W 



cos0-tj_sin0 = EcosO (100) 



-tj_ cos 



V± W 



EsinO (101) 



The above equation has two eigenvalues E± = — ¥± ± 



and clearly t± favors a state with a delocal- 



ized particle on the dimer. Thus in the ground state 



1 L ?> v 
cos = — — ,/H 



1 



V2 ] 



'1 - 



(102) 
(103) 



This state breaks the sublattice symmetry and there are 
clearly more particles on layer 2 than on layer 1. Sim- 
ilarly, if we deform ii) in analogous way, we will find a 
state with more particles on layer 1 than on layer 2. Both 
of these states are gapped. 

For infinitesimal i, we expect the energy of the broken 
symmetry state to be further lowered via second-order 



processes. This leads us to the conclusion that in the 
strong coupling limit, our Hamiltonian % has a ground 
state with broken inversion symmetry, i.e. the total num- 
ber of particles on the upper layer is different from the 
total number of particles on the lower layer. 



VI. SPIN-i CASE 

The symmetry based reduction of the number of cou- 
pling can be used for the spin-^ case as well. All the ar- 
guments presented in the section dealing with short range 
interactions follow through, but now the Fierz vector is 
18-dimensional, instead of 9. Specifically, each term in 
Eq. (|49[) . when multiplied by the appropriate coupling, 
gives rise to two terms as 



g X] (^X^) 2 -> gf^X^f + g^ a X j u a ^pf, 

(104) 

where the Pauli a corresponds to spin-i SU{2). The 
seemingly independent couplings in the two different 
channels, c and s, are still related to each other via a 
Fierz-like identity. 

In particular, we can use the SU(8) algebraic identity 



XnT,,,,, = ^Tr[SA a TA b ]A^A^ 



(105) 



where S and T are 8x8 matrices, and the 64 genera- 
tors A a can be obtained from the 16-SU(A) generators as 
{r a l,T a <8>cr z ,r a ® cr x ,r a ® a v }. This leads to 

_l Tr [SA a TA b ] (^( x ) A <V(y)) (^(2/)A<V(:#06) 

Again, the minus sign comes from tp and being anti- 
commuting (four component) Grassman fields. For con- 
tact terms x = y and the above equation constitutes a 
set of linear relations among the 18 symmetry allowed 
terms. 

If we now arrange the quartic terms into an 18- 
component vector V (Eq |Bl[) we can write the above con- 
straint as 



(107) 



TV = 



where the matrix displayed in Appendix B (|B2[) . has 
nine zero eigenvalues, and, as a result^! there are 9 inde- 
pendent couplings in the spin-i case. From here it can 
be shown that one can eliminate all the (spin-spin) g s x 
couplings in favor of the (charge-charge) g c x couplings 6 . 

Using the same technique as described for the spin- 
less case in the Appendix, we find the RG flow equations 
for the nine couplings in the spin- 1/2 case. These equa- 
tions, (| A10AA17]) . are shown explicitly at the end of the 
Appendix A. While full analysis of the Eqs. (IA10IIA17|) 
is beyond the scope of this paper, we have studied the 



tjj^Xj ® lip 


1 


2 


3 


4 


A A 


0,0,0,0,0,0,0,0,0 


1,-1,-8,-2,0,1,-1,2,0 


0,0,0,0,0,0,0,0,0 


1,-1,0,0,0,-1,1,0,-8 


A b 


1,-1,0,2,-8,1,-1,-2,0 


0,0,0,0,0,0,0,0,0 


0,0,0,0,0,0,0,0,0 


1,-1,0,0,0,-1,1,0,-8 


A c 


2,2,-4, -4,-4,2,-14, -4,8 


1,-1,0,2,-8,1,-1,-2,0 


2,2, -4,-16,4,-2,-2,0,0 


1,-1,0,0,0,-1,1,0,-8 


A D 


1,-1,-8,-2,0,1,-1,2,0 


2,-14,-4,4, -4,2,2,4,-8 


2,2,-4,-16,4,-2,-2,0,0 


1,-1,0,0,0,-1,1,0,-8 


Tp^Xj (g> ai/j 


1 


2 


3 


4 


A A 


0,0,0,0,0,0,0,0,0 


1,-1,0,-2,0,1,-1,2,0 


0,0,0,0,0,0,0,0,0 


1,-1,0,0,0,-1,1,0,0 


A b 


1,-1,0,2,0,1,-1,-2,0 


0,0,0,0,0,0,0,0,0 


0,0,0,0,0,0,0,0,0 


1,-1,0,0,0,-1,1,0,0 


Ac 


2,2,-4, -4, -4,2,2,-4,8 


1,-1,0,2,0,1,-1,-2,0 


2,2,-4,0,4,-2,-2,0,0 


1,-1,0,0,0,-1,1,0,0 


A D 


1,-1,0,-2,0,1,-1,2,0 


2,2,-4,4,-4,2,2,4,-8 


2,2,-4,0,4,-2,-2,0,0 


1,-1,0,0,0,-1,1,0,0 



TABLE II: The susceptibility coefficients (spin- 1/2 case) 
Axj , B Xj , C Xj , T) Xj , Sx 3 , Fx, ,Gx 3 , W-Xj , and X X] in Eq. {JSTJ 
for different particle-hole order parameters ip'Oiip. 



effect of three of these couplings in Ref^., starting with 
g^l, and generating and g^'. No other couplings 
are generated, assuming they vanish to begin with, in 
agreement with Eqs. (|A10IIA17|) . The equations pre- 
sented in the Appendix of this paper reduce to Eqs.(6- 
8) of Rcf. 3 provided wc set N = 4 there and identify 

.9i <-> 9aI' 92 <-> go] and 53 <-» 9a 2 - In this case : 

(c) . ... (c) 

for finite initial g A > and vanishing initial g\ 2 and 



gfo' 2 , the most dominant divergence appears in the ne- 
matic channel, which corresponds to one of the ratio sinks 
9*^d 2 / 9*aI = 7711 ~ —0.525. For certain combinations of 
the couplings, a different sink (9*$ 2 /9*aI = m 3 ~ 13.98, 
top fixed point in Fig. 3 in Ref<£) may be reached. Thus, 
an anomalous quantum Hall state may in principle be 
stabilized in weak coupling as well. 



Finally, we note in passing that the number (9) of in- VII. HUBBARD MODEL ON THE A-B 

dependent couplings in the spin-± case is in agreement STACKED HONEYCOMB BILAYER 

with Ref£, but disagreement with Ref A 



A. Susceptibilities 

Just as in the spinless case we can analyze the flow of 
various source terms A in order to determine the most 
dominant weak coupling ordering tendencies. Since in 
the spin- 1/2 case there are 9 independent coupling con- 
stants, we have 



In this section we use the above machinery to study the 
weak and strong coupling limits of the repulsive Hubbard 
model on the A-B stacked honeycomb bilayer. Just as 
before, we assume nearest neighbor hopping only, and 
the potential energy term can be written as 



din A x } 
din s 

+£x j gB 1 + Fxj9B 2 + Q Xj gc x + U X] g a +ix 3 gf3) ^~ 



2+(A Xj 9A 1 + B Xj g D2 + C X] g A2 + T> X] g 7 

> m 

(108) 



The coefficients A — I in 32 different particle-hole chan- 
nels are listed in Table HU The most dominant instabil- 
ity channel, Xj , yields the largest right-hand side of the 
above equation. 



v[H) = ^EE4r( R )Mft)«]i(RKi(R) 

3 = 1 R 

+ UJ2 &t it( R + ^)M R + ^) & u( R + S)b n (K + 5) 



R 



UJ2 b\ t (R - 5)b lt (R - ^)6 f u (R - £)& U (R - S) 

(109) 



R 



din A x 
dhis 
0.75 -, 



O.r, 



O.2.", 




FIG. 5: (Color online) Susceptibility vs In s (Eg |108|) for the 
Hubbard model with initial A m A U = 0.01, in the channels 



C} SJ - B\ a> (left to right). The fastest diver gence appears m 
the antiferromagnetic channel C\ ® S. Altogether 32 particle- 
hole channels have been analyzed (Table Hip : the channels 
not shown are either symmetry-related to the ones shown, or 
dlnA/dlns — 2 vanishes (or is negative). 



A. Weak coupling limit 



ergy effective Lagrangian can be written as 

U f ,o T / , + . - , \ 2 /,+ „ . , \ 2 



(111) 



\ (tfx* ® 1^) 



X=A,B,C,D 



This means that, of the 9 symmetry-allowed coupling 
constants, the only non-zero ones are g^, g^ 



with initial values 



and gfi', 



9<g(*=l) = ffg(*=l)= ' 



4A, 



flj c) (« = l) = 



U 
8A UC 



(112) 
(113) 



Next, we numerically solve the RG flow equations 
(|A10|) - (|A17|) with the above initial conditions (for 
itiU/(4:TtA uc ) = 0.01), and substitute the resulting s- 
dependent couplings into the susceptibility flow equation 
(jl08[) using the coefficients displayed in Table HU Com- 
parison of the resulting susceptibilities in 32 particle-hole 
channels shows that the most dominant divergence ap- 
pears for O = C\ ® 3. Physically, this corresponds to an 
anti-ferromagnetic state, with anti-aligned spins on the 
sites b\ and 62. 



Projecting the Hubbard interaction onto the low en- 
ergy modes we find 



where 
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(110) 



The 8-component Fermi fields ip are understood to be at 
space- (imaginary) time point r, r. 

Using the notation established in the previous sections, 
we note that the above interaction terms in the low en- 



B. Strong coupling 

It is well known that— in the strong coupling limit 
[//t > 1 the Hubbard model with one particle per site 
is equivalent to the spin- 1/2 Heisenberg model with an- 
tiferromagnetic coupling J ~ t 2 /U. If we set t± = 0, 
the two honeycomb layers decouple and at strong cou- 
pling each layer orders anti-ferromagnetically 2 ^ with a 
sublattice magnetization that is free to point along any 
direction on each layer. Once t± is finite, the sublat- 
tice magnetizations on the two different layers lock into 
relative anti- ferromagnetic arrangement. 

We thus find that the half-filled Hubbard model on the 
A-B stacked honeycomb bilayer orders antiferromagneti- 
cally in both the weak and strong coupling limits. 



VIII. CONCLUSIONS 

We have studied the effect of short range interactions 
on fermions moving on the A-B stacked bilayer. In or- 
der to access the "strong coupling phases"— from weak 
coupling, we have fine-tuned the spectrum of the non- 
interacting Hamiltonian to achieve parabolic degeneracy, 
with the ensuing logarithmically-divergent susceptibili- 
ties that appear in several channels. We have found that, 
in the spinless fermion case, the typical dominant order- 
ing tendency opens a spectral gap, although the nature of 
the resulting insulating state may be dramatically differ- 
ent. For example, the weak coupling limit of the spinless 



t — V model, with nearest neighbor hopping and near- 
est neighbor repulsion, leads to an inversion symmetry 
breaking (trivial) insulating phase, while the (right) sink 
in the RG flow diagram shown in Fig. ^ corresponds to a 
spontaneously time-reversal symmetry breaking anoma- 
lous quantum Hall phase with a xy = ±2e 2 /h. Under 
certain conditions, the dominant instability may appear 
in the particle-particle channel as well. In addition to 
the generic instabilities of the spinless model, fine-tuning 
of the initial couplings may lead to a flow towards the 
non- interacting fixed point. While such behavior is non- 
generic, it is interesting that there is an entire surface in 
the ratio space (FigEJ) which gives rise to such a flow for 
positive initial gA 2 ■ We have also studied the strong cou- 
pling limit of the t — V model. In this case, the ground 
state wavefunction can be shown to explicitly display in- 
version symmetry breaking. Since the same type of order 
is found in the asymptotic limits of strong and weak cou- 
pling, it is reasonable to assume that for this specific 
model, such ordering happens for any (repulsive) cou- 
pling strength. 

In the spin-1/2 case, we find 9 independent, symmetry- 
allowed couplings, and their RG flow equations. While 
these equations have not been studied in their entirety, 
they reduce to the ones presented before in Ref£, in 
which case an analysis similar to the one presented for 
the spinless case here, leads to either a nematic phase or 
an anomalous quantum Hall phase with a xy — ±4e 2 /h, 
where the extra factor of 2 compared to the spinless case 
is due to trivial spin degeneracy. 

Moreover, these equations (|A10|) - (IA17|) are solved nu- 
merically for the spin-1/2 Hubbard model at half fill- 
ing. The initial values of the effective couplings are such 
that the most dominant particle-hole instability appears 
in the anti-ferromagnetic channel. This dominance has 
been established by comparing susceptibilities toward 32 
different ordering tendencies. Since the same instability 
appears in the strong coupling limit, it is reasonable to 
conclude that the antiferromagnetic order sets in for any 
U > 0. 

When trigonal warping is taken into account, the loga- 
rithmic infra-red divergences are cut off 3- ^ by the energy 
scale corresponding to the deviation from the parabolic 
spectrum. This means that the non-interacting system 
is stable towards infinitesimal coupling, i.e. there are no 
true weak coupling instabilities. Instead, the interaction 



strength must be increased beyond a critical value, which 
may be quite difficult to obtain accurately. However, it is 
worth noting that no such fine-tuning is necessary for the 
models with parabolic touching studied in Ref. ? , in which 
case there are true weak coupling instabilities. Unfortu- 
nately, one cannot just immediately translate the results 
regarding the dominance and the nature of the weak cou- 
pling broken symmetry states found here for the honey- 
comb bilayer, because 1) the location of the degenerate 
points in the Brillouin zone is (qualitatively) different and 
2) the lattice symmetry will, in general, allow different 
contact terms than those found here. 

Acknowledgments 

This work was supported in part by NSF CAREER 
award grant No. DMR-0955561. 



Appendix A: Details of the RG derivation 

For general coupling constants gsr, expanding in pow- 
ers of g, gives the cumulant expansion 

U-kasT S^S^t^k _ e -ig ST (/ 1 i/.ts^ t rv(i)) 
xexp \ 9ST9UV [ (((^SrJj^Tip(l)) (</> f W f V^(2))) 

L 8 Jl,2 

- ((^Sip^Tip(l)))((^Utp^Vtlj{2))))] 

where the average (. . .) is with respect to the gaussian 
weighting factor. We have used a short-hand 1, 2 for the 
modes at the space- (imaginary) time 7^2, 1*1,2 and each 
ip = V ; > + ip<- We integrate over the fast modes ■0> 
whose wavenumbers A/s < k < A. The non-interacting 
Green's function is 

G k M = (-^iS-rf,)' 1 = *" + S ' 2 dk 2 (Al) 
and just as before d? = — -, dl — ^ Ji , and, in the 

J K 2m ' K 2m 5 ' 

spinless case, T, x = 72 and S y = 71. 

Using the identities ([58]l . we can evaluate the needed 
diagrams. All possible contractions correspond to the 
diagrams in the Figure ©. 



For the first diagram we find the following terms 



AS { e ff A} = ±££ 9S9U I I V f (l)^(l)Tr [SG(1 - 2)UG(2 - 1)] ^(2)U^(2), 



see ueg 



(A2) 



where, in the spinless case, Q = {A\, A%, Di, D2, C3, D3} and the corresponding couplings, in order of appearance 
of S in Q, are {gAn 9a 2 > 9a 2 t 9d 2 i 9-yi 9-y}- Using the gradient expansion to determine the RG fate of the marginal 
couplings, we find that 



AS 



(RPA) 



\Y,Y.9s9u I ^ f (l)^(l)Tr 
segueg Jl 



-SU + -SyiUji + -S-y 2 U^f2 



i>\l)U^{\)— Ins. (A3) 

47T 



Performing the traces gives 



AS<* PA) =-j (2g% [{^A^f + (^D^) 2 ] + 4<4 (^ D 2 ^f + Ag 2 [{^C^f + D^f] ) £ Ins. (A4) 

For the second and third (vertex) diagrams in Figure [2] we have the following terms: 

AS %] = - E E 9S9U [ [ 4>Hl)SG(l - 2)UG(2 - 1)^(1)^(2)^(2). (A5) 
sesc/es 

Performing the gradient expansion gives 



AS S = [ ^(1)5 (-f/+^7i^7i + ^72C/72) ^(1)^(1)^(1)^ Ins. 

see [/eg ' /l v 



(A6) 



Performing the requisite sums and matrix algebra gives 

AS*}] = -^"^(-g^ +.gc 2 +2.g 7 ) [(^^^(V^i^) 2 ] -2 ffD2 (5A 1 -2 5 A 2 +5D a +2 ff7 )(V t r>2i/') 2 )£lns 
- J (-2.g 7 ( 5Al - 2g A2 +g D2 ) [{^C^f + (V^) 2 ]) ^ Ins. (A7) 

For the fourth and fifth diagrams in Figure [2] we find the following terms: 

AS e L fl = -7 E ^ / ((V' t (i)[5,^V'(i)) 2 + ^E(^ t ( 1 )( 5 ^ C7 +^)^ 1 )) 2 ) £ lns - ( A8 ) 

S,U£Q Jl V o=l / 

The corresponding change in the effective action is 

= - (2g Al g A2 (^A^) 2 + ^(51+4^-4^^+2^+2^) [(V^) 2 + D^) 2 ]^ 
(2g A2 (g D2 - g A2 ){^D 2 ^) 2 + 2g A2 g y [{^ C^) 2 + D^f] - 2g 2 (^B 2 ^) 2 ) 
((g D2 - 2g A2 )g 1 [{^ A^f + {^B^f + {^G^f + (^D^) 2 ]) . (A9) 



In order to find the renormalization of the coupling constants, the last two terms must be rewritten using the Fierz 
identities (|52|) and (f56|) . Adding the terms from AS^f A \ AS^JJ, and AS^j, rescaling the fields and the integration 
measure, and comparing to the starting action (|57p we find the RG equations (|60[) for the four coupling constants in 
the spinless case. 

The above equations (IA2[) . (|A5I) and (|A8|) can also be used to derive the flow equations of the 9 coupling constants 
in the case of spin-1/2 fcrmions. In this case, we have Q = {Ai ® 1, A% (E> 1, Z?i ® 1,-Bi S3 1, C2 ® 1, £?2 <8> 1, G\ ® 1, Z?2 (8> 
1, A3® 1, i?3<8) 1, S4® 1, Aa® 1, C4® 1, £>4<8> 1, (73® 1, £>3<g> 1} and the corresponding couplings, in order of appearance 

01 o m y, are \g Al , g A2 , g A2 , g Bl , g Bl , g B2 , g Cl , g D2 , g a ,g a ,yp ,3^ ,33 i9p t9i 7S7 /• 



The resulting RG flow equations for spin-1/2 fermions are 
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dins 



<ilns 
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til +*9%9%-Wi 
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47T 



2 (V + (j 

fag (*2 - %g - 4 9 g + 9 g - 9 g + flg - 2. 9 i c ) + 2gM 
2 (>gsg - $(2*2 - ffg) - 2<£>(2<#> - #>))) £ 

-4 (Ag + *2*2 - ^ C)2 - 2 # )2 + - *? ,a ) £ 

4 ( 9 g („« - 3| ,2 ~ 2 5 g + 9 g _ 35 g + fl g _ 2g ( C ) + 4 5 M _ 2^) 
. 9 g(2^- 5 g)-2^( 5 g-^))£ 
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47T 
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(c) _ „(c) 

a 2 ys 2 



(A12) 

(A13) 

(A14) 
(A15) 



(A16) 



-4 ( 5 i c) (,g - ff g) + 4 C) (*2 - 2.2 + - ( 5 g - 3 ff g + 2sg - fl g + fl g - % 
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Appendix B: Fierz reduction for spin-| case 

If we arrange the quartic contact terms into an 18-dimensional vector 
V = {tylAM 2 , (^A 2 V«) 2 + {^ a D^ a )\ (^ a B^ a ) 2 + {^ a C 2 ^ a f, tyiBM 2 , (^ a C^ a ) 2 , 

(ipiD 2 a a piljp) 2 , (iplA 3 a a i3ipi3) 2 + {^ a B z a a pi)p) 2 , ^ (^i x i^ai3^p) 2 , (i>iC 3 a a piJ)p) 2 + (ipiD 3 a a pipp) 2 , > 

X=A,B,C,D J 

(Bl) 

then the Fierz identities (| 106[) can be used to relate different components of V via the linear constraint TV = 0. 

In practice, to obtain the Fierz matrix T we arrange the 64 SU(8) generators A a in the order {X <g> 1,X (g> 
CT z ,X(Kicr ;r ,X(g)CT 2 '}, where X = {A x , A 2 , D 1 , B 1 , C 2 , B 2 , C\, D 2 , A3, B 3 , B 4) A 4 , C 4 , D 4 , C 3 , £> 3 }. Straightforwardly, a 
64-component analog of the Fierz vector, obtained by using our ordered set of SU(S) generators, will be denoted by 
V. Next, we notice that only the (diagonal) terms with S = T enter our Cmt, and that Tr[SA a 5A b ] ~ S a b. Therefore, 
for A a arranged as described above, we numerically generate a 64x64 matrix $ a b = — Tr[A a A 6 A a A 6 ]/64. We then 



construct an auxiliary 64x64 matrix 
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It is easily seen that components 1 — 9 and 17 — 25 of MV correspond to our Fierz vector V. Moreover, the blocks 
1 — 9 and 17 — 25 of the matrix M&M -1 do not couple to the rest of the components, and, when subtracted from an 
18-dimensional unit matrix, correspond to the sought Fierz matrix, 
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in order to write 




MV. 



One can check that the above matrix has 9 zero eigenvalues, which implies 9 independent couplings and 9 constraints. 
In addition, one can solve the linear system TV = and eliminate all terms of the form (ip^Xjaapipp) 2 in favor of 

($Wa) 2 . ! 
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